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Modeling of thermal transport in practical nanostructures requires making trade-offs between the 
size of the system and the completeness of the model. We study quantum heat transfer in a self- 
consistent thermal bath setup consisting of two lead regions connected by a center region. Atoms 
both in the leads and in the center region are coupled to quantum Langevin heat baths that mimic 
the damping and dephasing of phonon waves by anharmonic scattering. Increasing the strength of 
the coupling reduces the mean free path of phonons and gradually shifts phonon transport from 
ballistic regime to diffusive regime. In the center region, the bath temperatures are determined 
self-consistently from the requirement of zero net energy exchange between the local heat bath and 
each atom. By solving the stochastic equations of motion in frequency space and averaging over 
noise, we derive the formula for thermal current, which contains the Caroli formula for phonon 
transmission function and reduces to the Landauer-Biittiker formula in the limit of vanishing bath 
coupling. We prove that the bath temperatures measure local kinetic energy and can, therefore, 
be interpreted as true atomic temperatures. In a setup where phonon reflections are eliminated, 
Boltzmann transport equation under gray approximation with full phonon dispersion is shown to 
be equivalent to the self-consistent heat bath model. We also study thermal transport through 
two-dimensional constrictions in square lattice and graphene and discuss the differences between 
the exact solution and linear approximations. 
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I. INTRODUCTION 

Recent theoretical and experimental studies of ther- 
mal properties of materials have demonstrated many ex- 
otic phononic phenomena such as ballistic and anoma- 
lous transport [iHl , conductance quantization and 
phonon tunneling 0| . These discoveries suggest that 
the ability to manipulate heat flow at microscopic level 
and to better understand phonon transfer in nanoscale 
may lead to important technological breakthroughs rang- 
ing from new materials for thermoelectric conversion 
Q to improved thermal management in future elec- 
tronics [loj and even information processing by phonons 

[IMl. 

Modeling of thermal transport in practical nanostruc- 
tures typically requires making trade-offs between the 
size of the system and the completeness of the the model. 
Consequently, the commonly used models such as Boltz- 
mann transport equation (BTE) [l^ . molecular dynam- 
ics (MD), Landauer-Biittiker (LB) formalism [isl . [l6| 
for phonon transfer 0, Il7| and full non-equilibrium 
Green's function (NEGF) method [H, [H] each have dis- 
tinct strengths and weaknesses. For instance, BTE for 
phonons is a powerful method that is applicable even for 
macroscopic systems, but it does not apply well to micro- 
scopic systems where wave effects such as diffraction are 
important. MD can be applied to phonon transport in 
microscopic systems and accounts, e.g., for wave effects 
and phonon-phonon scattering due to anharmonicity of 
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the interatomic potential, but it becomes computation- 
ally heavy for large systems and cannot strictly account 
for quantum statistics. The LB and NEGF models can 
fully account for the quantum statistics, but LB assumes 
ballistic phonon transfer and NEGF is computationally 
extremely demanding and therefore limited so far to very 
small systems (20| . 

As a consequence of the above limitations, none of the 
above models are well suited for modeling phonon trans- 
fer in typical nanostructures consisting of a relatively 
large number of atoms. A very interesting compromise 
between system size and model completeness is provided 
by the self-consistent thermal bath (SCTB) model sug- 
gested by Bolsterli, Rich and Visscher (BRV) [2l|. In 
the SCTB model, the phonon scattering is mimicked by 
coupling the atoms to local heat reservoirs whose tem- 
peratures are determined from the condition that, in the 
steady state, there is no net energy transfer between an 
atom and the corresponding local heat reservoir. The 
concept was first used to show that for a classical system 
with bath temperatures equal to the local kinetic tem- 
peratures the thermal conductivity of a harmonic one- 
dimensional chain was rendered finite by the bath cou- 
plings. More recently, the self-consistent thermal bath 
model has been applied to studying local thermal equili- 
bration (2^ , quantum effects in non-ballistic heat transfer 
[23l - l25t and the necessary ingredients of thermal rectifi- 
cation [H^. 

In this paper, we extend the SCTB models beyond 
one-dimensional chains and study the heat transfer and 
the use of SCTB models in describing quantum ther- 
mal transfer in one-dimensional and two-dimensional 
structures that exhibit geometric as well as phonon- 
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phonon scattering. We compare the predictions of SCTB 
with Landauer-Biittiker (LB) formahsm and Bohzmann 
transport equation (BTE). In contrast to LB formahsm, 
phonon transport in the SCTB model is not purely bal- 
listic due to the interaction with the local heat baths, 
but we show that the SCTB model reduces to the con- 
ventional LB model in the limit of vanishing bath cou- 
pling. In a setup where wave effects can be neglected, 
SCTB is shown to be equivalent to BTE under gray ap- 
proximation. We also demonstrate how the local bath 
temperatures are intuitively related to the local energy 
densities. 

The paper is organized as follows. In Sec. [ITl we 
present the computational setup and derive the for- 
mula for heat currents flowing to the leads and to self- 
consistent heat baths. This is achieved by first solving 
the Heisenberg-Langevin equations of motion in Sec. Ill Al 
and then specifying the statistical properties of noise 
terms in Sec. IIIBI Sections III CI and IIIDI are devoted 
to calculating the average heat flow to the baths and 
presenting physical interpretation for the self-consistent 
bath temperatures. For comparison purposes, we also 
solve the Boltzmann transport equation under gray ap- 
proximation for the one-dimensional chain in Sec. Ill El 

In Sec. mil we discuss how to solve the self-consistent 
equations either by iterative means or by linearization. 
We also present a physically intuitive method of solving 
the equations, which can be interpreted as describing the 
transient behavior of the system. As an application of 
the formalism, we study in Sec. IIV Al heat transfer in a 
one-dimensional chain coupled to semi-infinite chains, in 
the so-called Rubin-Greer setup [1^, and highlight the 
connection to Boltzmann transport equation. We then 
study thermal transport through two-dimensional con- 
strictions both for square lattice and the more practical 
case of graphene in Sees. lIVBI and llVCI The methods for 
solving the self-consistent non-linear equations are com- 
pared in Sec. IIVDI Conclusions are finally given in Sec. 

El 



II. THEORY 

In the theory and results of this paper, we mainly 
focus on a system that essentially consists of the left 
lead region, center region and right lead region as shown 
schematically in Fig. [1] All atoms within the leads are 
coupled to local Langevin heat baths set to prescribed 
values Tl and Tr. The atoms in the center region are 
coupled to local heat baths whose temperatures are de- 
termined self-consistently from the requirement of lo- 
cal current conservation. The coupling to the Langevin 
heat baths effectively mimics thermalizing events such as 
phonon-phonon scattering. It is important to stress that 
in contrast to Landauer-Biittiker model, phonon trans- 
port is not assumed to be ballistic either in the leads or 
in the center region. 

In the following, we first solve the equations of motion 
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FIG. 1. (Color online) (a) A schematic illustration of the 
system under study. The structure is divided into the left 
lead, the center region and the right lead. All atoms are cou- 
pled to spatially uncorrelated quantum Langevin heat baths, 
which are shown explicitly for one cross-section in (b). In 
the left and right lead, the temperatures of the local heat 
baths have prescribed values Tl and Tr, respectively. In the 
center region, on the other hand, temperature varies between 
Tl and Tr and the bath temperatures are determined self- 
consistently using the requirement that the average thermal 
current to each bath vanishes. The leads can contain an in- 
finite number of atoms, but the center region is finite. Two- 
dimensional square lattice with nearest neighbor interactions 
is shown for illustrative purposes, but the basic principle can 
be applied to any geometry. 



in Sec. Ill Al then specify the statistical properties of the 
noise in Sec. Ill Bl and finally derive the formula for heat 
currents in Sec. Ill CI 



A. Equations of motion 

The time evolution of atoms in the setup of Fig. [T] 
consists of two parts. The first part is deterministic and 
is specified by the system Hamiltonian Tl and Heisen- 
berg equations of motion. The second part consists of a 
stochastic force and friction due to the interaction with 
the local heat bath and cannot be directly derived from 
a Hamiltonian [soj . 

The Hamiltonian time evolution of the atomic displace- 
ment uf of atom i along direction a £ {x, y, z} and 
corresponding conjugate momentum pf is determined by 
the Hamiltonian H and the Heisenberg equations of mo- 
tion iif = {i/h)[H,uf] and pf = {i/h)[H,pf]. Here 
[A, B] = AB — BA denotes the commutator and the 
atomic displacement uf = qf — qf'^ is defined as the vari- 
ation of position qf from the equilibrium position qf^. 
The Hamiltonian of the system is, in the harmonic ap- 
proximation, 



n 



2 ^ 



Uj K/U/ 



|EE"FV,.u,, (I) 



where index / S {C, L, R} labels the region: C stands for 
center region, and L and R for the left and right leads. 
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respectively. The displacement and momentum vectors 
u/ and Pi contain the displacements and momenta of all 
particles in region / and we assume the masses m of all 
atoms to be equal for notational simplicity. The spring 
constant matrix K/ and the coupling matrices V/j are 
the block components of the full spring constant matrix 
K divided into blocks as 



K = VcL Kc VcR 
\ Vhc Kh 



(2) 



where we assumed that the leads do not interact, so 
Vlr = = 0- The elements of K are obtained from 
the second derivative of the interatomic interaction en- 
ergy V as (sH 



(3) 



The equilibrium positions are defined by the condition of 
zero force 



dV 



(4) 



which must be satisfied for all atoms i and components 
a. 

The Heisenberg equations of motion that follow from 
the quadratic Hamiltonian ([T]) coincide with the classi- 
cal equations of motion. Accompanied with the non- 
Hamiltonian time-evolution arising from the interaction 
with the heat bath, the equations of motion become 



-Kiui 



J2 ^i-"'-' 



The last two terms are Langevin friction and noise terms 
that turn the Heisenber g eq uation of motion into a quan- 
tum Langevin equation |24l. Isol. Is^l . The stochastic force 
^7 is a vector whose I'th component is the fluctuating 
force at site i due to the interaction with the local heat 
bath. The statistical properties of the Langevin terms 
are discussed in the next section. 

Focusing on the steady-state behavior enables solv- 
ing Eq. ^ by Fourier transformation defined, as usual, 
by /(w) ~ J dte'^'^* f (t) with the corresponding inverse 
transformation /(t) = J(daj/27r)e^'"*/(w). Equation (O 
transforms into 



mu!' 



■ui = -K/U/ - ^ V/jUj + imjjujui + ^j. (6) 



Rearrangement of Eq. (O gives 



u/(w) = g/(w) ^Vijuj{u}) - ii{uj) , (7) 
j=ii J 
where the uncoupled Green's function is defined as 



g/(a;) = [mw^ — K/ -|- im^iuj] ^ . (8) 



The uncoupled Green's function includes damping self- 
energy irwyju) due to coupling to the heat baths. 

Substituting Eq. (0 for / = L, i? to dH) for / = C 



m7/U/+^/. (5) gives 



^uc(w) = -Kcuc(^^) - ^ Vc/g/(w) V/cUc(w) - |/(a;) -t- im7ca;uc(w) -f |c(w) (9) 

I=L,R 

= -KcUc{^)- [J:i{oj)uc{oj) - fiijoj)] + imjc(^uc{oj) + $c(^)- (10) 



I=L,R 



In the second line, we defined the lead self-energies 

J:j{uj) = Vc/g/(c^)Vjc (11) 
and the lead-coupled Langevin noise terms 

fii{u:)=VciSiiL^)ii{^)- (12) 
The solution to Eq. ^ is 



I=L,R 



(13) 



where the full Green's function of the center region is 



G(^) 



— K, 



c + irti'icui - ^ S/(a;) 



I=L,R 



(14) 

Equations (fT2|) and (fT3| state that thermal fluctua- 
tions in the leads can propagate to the center region as 
described by the Green's function g/(w) and coupling 
matrix V^/ and thereby introduce additional noise terms 
fji^ and fiji in the center region. The self-energies T,]^ and 
T,R appearing in the Green's function (fT^ describe the 
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energy shift and broadening of the phonon energy levels 
in the center region due to the leaking of phonons into 
the leads. The self-energies of the semi-infinite leads can 
be determined by using, e.g., the recursive decimation 
routine by Lopez and Sancho (33| . 

By integrating out the leads, we have effectively re- 
placed the lead coordinates by the noise terms 77/ and 
the accompanying self-energies S/. In the next section, 
we derive the fluctuation-dissipation relation connecting 
the statistical properties of fij to Im[E7(w)]. 

B. Noise power spectra 

The Langevin noise operators ^/ appearing in Eq. ([S]) 
act as stochastic sources of thermal fluctuations due to 
coupling to the local heat bath [U, [s^, H^]. They are 
Gaussian random variables with zero mean and covari- 
ance related to bath temperature. To calculate the sta- 
tistical averages of observables such as heat current in 
the center region, we need the covariances of both the 
bare Langevin noises ^-nd the lead noise contributions 
77/ = Vcigiti- Since the bath temperatures in the center 
region depend on position i, we handle the noise terms 
originating from the center region and leads separately. 
In the following, we assume for notational convenience 
that each atom only has a single degree of freedom corre- 
sponding to displacement along, say, a;-direction. In the 
general case, local bath at site i will be coupled to dis- 
placements {uf ^u}^ ,ul) in different coordinate directions 
with a single temperature T!j, making the notation a bit 
more cumbersome but analogous. 

The covariance of the noises produced by the local heat 
baths at sites i and j in the center region is, for h = ks ^ 

1, MM 

where the coupling function for Ohmic friction is 
Tij{uj) = 2m'yc'-^Sij. This corresponds to the memoryless 
friction assumed in Eq. ((S]), but more general couplings 
could be straightforwardly included as well. The friction 
parameter 7c determines the strength of coupling to the 
heat baths and can be interpreted as phonon decay rate 
[s^ l- The corresponding scattering time is tc = ^c^, 
which is independent of frequency for Ohmic baths. 

The term in braces, where the bath temperature 
Ti appears in the Bose-Einstein function fsiujjTi) = 
[exp{uj/Ti) — can be written in the more trans- 

parent form fB{uJ,Ti) + 1/2 -|- 1/2, where the first term 
is the thermal phonon occupation number, the second 
term comes from zero-point fluctuations and the last 
term reflects the non-commutative quantum nature of 
the Langevin operators [IHl- One can write the term as 
the sum of odd and even functions in a; as (w, T,) -|- 1 = 
coth(w/2ri)/2 + l/2, and it turns out that the additional 
factor of 1/2 cancels in all integrals over frequency after 
proper symmetrization. 



For the noise operators in the leads, the temperatures 
of the baths have prescribed values Tl and Tr, which 
do not depend on position. Therefore, we can write the 
covariance directly in matrix form as (/ £ {£, R}) 

{il{iu)ij{w'f) = 27T6{UJ+Uj')r'{uj) [fB{LO,Tl) + 1] , 

(16) 

which is useful in the following calculations. Here, the 
coupling function matrix is a diagonal matrix with 
elements r^^ (a;) — 2mjjuj5ij. 

Using iiiiuj))) = and Eqs. (HU and (HH), we see 
that the noise terms fji originating from the leads satisfy 
(77/ (w)) = and 

{fji{io)fii{Lu'f) = 27r(5(w + Lu')liiLj), (17) 
with the power spectrum 

hiiu) = Vcigi{uj)f'{L,)gi{^uj)Yic [fB{oJ,Ti) + 1] , 

(18) 

where we noted that the Green's function gi{uj) is sym- 
metric, since the spring matrix K/ is symmetric. A 
straightforward calculation shows that 

SO Eq. (fT8)) becomes 

I/(w) = lYci (giiuj) - gi{cjy)Vic [fB{uJ,Ti) + 1] 

(20) 

= -2Vc/Im [giiu)] Vic [/b(^, Ti) + 1] (21) 
= -2Im[I]/(w)] [fB{io, Ti) + 1] , (22) 

where we used the definition (|lip . Defining the lead cou- 
pling function 

r^{uj) = -2Im[I]7(w)], (23) 

we see that the power spectrum of the noise caused by 
the leads can be written as 

(fji{io)f)iiuj')'^) = 2TrS{uj + uj')r' {io) \fB{oJ,Ti) + 1] . 

(24) 

Equation ([24| is analogous to Eqs. (fT5|) and (fT6|) ex- 
cept for the form of the coupling matrix r^(a;), now de- 
fined using the self-energy of the lead as shown in Eq. 
([23]) . Equation ([24|) is one of the main results of this 
paper, showing that an atomic reservoir (lead) coupled 
to local heat baths at prescribed temperature can be 
represented by noise and dissipation terms related by a 
fluctuation-dissipation relation. The generalization com- 
pared to Ford-Kac-Mazur formalism [24, 36, is that 
we do not assume heat transfer in the reservoir to be fully 
ballistic. 

Solution (fT5|) combined with the noise correlations (fT5)) 
and (|24p allows us to calculate the thermal averages of 
all observables of interest. 
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C. Heat flow to baths 

In the Heisenberg-Langevin equation of motion (fTO| . 
the friction and stochastic force terms induce energy ex- 
change with the heat bath. The energy exchange rate can 
be calculated from the time derivative of local energy. A 
natural definition for the local Hamiltonian of atom i in 
the center region is 

^' = ^j + IJ2''^^^^''^- (25) 
j 

In Eq. p5p and from now on, we drop the index C de- 
scribing the center region, since the lead coordinates do 
not appear anymore. Using the equation of motion ([S]), 
the symmetrized time derivative taking into account the 
non-commutativity of Ui and pi can be calculated to be 

^» = ^ ml + i^-^'J' ({"*'"J^ + {ui,Uj}) (26) 
j 

j 

Here {A, B} = AB + BA is the anti-commutator and 
for simplicity, we have assumed that the particle is not 
at the boundary so that it is not directly coupled to the 
leads. The term inside the sum in Eq. (|27p is the heat 
current flowing from site i to site j and the second term 
in parentheses is the heat current 

Qi = mjiif - i [ui£,i + £,iUi] (28) 

flowing to the local heat bath at site i. 

As shown in the appendix, the statistical average of the 
heat current is formed as a sum over the contributions of 
the left ( J = L) and right (J = R) leads and each local 
heat bath ( J s {1, 2, . . . , Nc}) as 

(Q,)= / |^27m^2^[G(c.)r^(c.)G(-c.)^]^^ 

X [fBiio^Tj) - fB{Lo,T,)]. (29) 

Here the sum over bath index J £ {L, i?, 1, 2, . . . , Nc} 
separately accounts for the contribution of each individ- 
ually treated heat bath to the thermal balance at site i 
as detailed in the appendix. The coupling matrix F"^ is 
defined by Eq. (|23| for the lead heat baths (J = L or 
J = R). For the local heat baths ( J e {1,2,..., Nc}), 
the only non-zero element of the coupling matrix F"^ is 
rjj = 2^muj. The term [GF'^G''^]i,; describes the ther- 
mal coupling between bath J and site i. In the following, 
we refer to the bath at site i simply as bath i. 

Using the definition of F* for local heat baths, we can 
write Eq. for the heat fiow to bath i in the general 



form 

/•OO T 

(Q.) = J ^^E'^^(^) ifBi^'Tj) - fB{uj,m , 

(30) 

where the transmission function between baths i and J 
is 

- Tr [F'(c^)G(c^)F'^(l.)G(-c^)] . (31) 

Equation ([50]) is also valid for the currents flowing to the 
leads, i.e. for the substitution z — > L or i — > i?, and the 
derivation proceeds analogously. In this case, one should 
use the equation of motion pO|) in the pi term of ((27|) to 
calculate the heat in-flow to center region by the noise 
term ?)/ and out-flow by the force term I]/u(a;). 

Equation pop is the multiprobe Landauer-Biittiker for- 
mula [l^ for thermal transfer between several heat baths. 
Equation (|3ip is the Caroli formula for phonon trans- 
mission function, flrst derived by Mingo and Yang (l7| 
from the mode picture and by Yamamoto and Watanabe 
(3^ using Keldysh formalism. We have rederived the for- 
mula using local Langevin heat baths and thereby also 
included dissipative effects in the leads. 

We point out that although the average heat current 
{Qi) vanishes in the self-consistent temperature conflgu- 
ration for all local heat baths in the center region, the 
spectral heat current 

mu;)) = 27™^.^^ [G(c.)F'^(c.)G(-c.)^]^. 
,7 

x[fB{Lu,Tj)-fB{L^,T,)] (32) 

to a local heat bath is generally non-zero. For example, a 
bath may have a net in-fiow of high-energy phonons, but 
then there must be a corresponding net out-flow of low- 
energy phonons. These non-zero spectral currents lead to 
the redistribution of phonon energies inside the structure, 
similarly to the full non-equilibrium Green's function for- 
mahsm where generally {Ql{uj)) ^ —{Qr{u)) jT8| . 

In the limit of vanishing bath couplings, 7 — >■ 0^, 77 -4- 
0"*", the lead and center region Green's functions ([8]) and 
region reduce to their ballistic counterparts and the 
only non-zero transmission function is Tlr{uj). Equation 
([30| reduces to 

poo J 

{Qn) = / -^^^.Tr [T^{u)G{u)T^{uj)G{-u)] 
Jo 27r 

x[!b{lo,Tl)- iBiLo^Tn)] (33) 

for the current flowing to the right reservoir, which is 
the two-probe Landauer-Biittiker formula for ballistic 
phonon transfer jT^ . 

D. Physical interpretation of the bath 
temperatures 

For the classical self-consistent thermal bath models, 
the requirement of zero net energy exchange with the 
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local baths was enforced by requiring the bath tempera- 
tures to be equal to the local kinetic temperatures (2H . 
The present model allows finding a fully quantum inter- 
pretation for the self-consistent bath temperatures. To 
this end, we first note that the heat current fiowing to 
a self-consistent reservoir at site i can be written in the 
form 



duj 
2^ 



uiDi{lo) 



fB{io,T,) 



(34) 

where the local kinetic energy is ef™ = mu|/2 and the 
local density of states (LDOS) is defined as Di{uj) = 
—4:Ujmlm[Gii{uj)] [40'|. This form results from noting that 
the first term of Eq. (^5]) is jmuf ~ 276^^™ , and the sec- 
ond term follows from Eq. (jA.10[l by using the definition 
of LDOS and dropping the odd term that cancels out in 
the integration. The self-consistency criterion (Qi) = 
then reduces to the requirement 



duj 
2^ 



ujDi{uj) 



/b(^,TO + 



(35) 



The left-hand side of Eq. ([55]) can be interpreted as the 
total energy at site i consisting of the kinetic and elastic 
energies, which are equal in a statistical-mechanical sys- 
tem according to the virial theorem [4l|. Virial theorem 
is, of course, rigorously valid only for the total kinetic and 
interaction energies at thermal equilibrium. The right- 
hand side is the total vibrational energy of an oscillator 
at temperature Ti. Equation (|35p gives a very natural 
interpretation to the self-consistent bath temperature Ti 
as a measure of energy located at site i. 

Note that for Ohmic baths, the integrals in Eqs. ([34|) 
and ([35|) actually diverge, because the density of states 
scales as D{u!) ^ — wlm[l/(w^ + i'^muj)] ^ uj^"^ for 
a; — >■ 00, resulting in a logarithmic divergence in the zero- 
point term. The divergence is, however, cancelled by an 
identical term in (ef™), making Ti well-defined. 

In the classical limit, Eq. ([55)1 reduces to 



27r 



(36) 



where we used the sum rule {duj/2Tr)Di{oj) = 1. This 
sum rule has been proven for the electronic case [i^ and 
the proof for phonons is analogous. Equation ([36)1 can 
be interpreted as the local equipartition theorem analo- 
gous to the statistical mechanical equipartition theorem 
(efcin) = where Nf is the number of degrees of 

freedom in the system. Relation (|36[) is routinely used 
as the definition of local temperature in classical molec- 
ular dynamics simulations [43j. Equation ([55)1 suggests a 
similar definition for quantum systems. 



E. Solution of the Boltzmann transport equation in 
ID chain 



In a sense, self-consistent thermal bath model (SCTB) 
can be thought of as the fully wave enabled extension of 



the gray approximation [ij, |4J| to the Boltzmann trans- 
port equation (BTE). Therefore it is instructive to com- 
pare the results obtained from BTE and SCTB under 
conditions where the wave-effects are negligible and there 
are no refiections between the chain and the reservoirs. 
For the simple one-dimensional string of length L, BTE 
in the continuum approximation reads 



dn+{x,uj) 7i+(a;,w) - no(a;,cj) 
v(u}) = (37a) 

ox T 

^ dn^{x,uj) _ n^{x,uj) - no{x,Uj) 
-v[u;) - ^ , i37bj 



where ?i-|-(a;,w) and n-{x,uj) are the distribution func- 
tions for states with positive and negative group veloci- 
ties, respectively. The thermal boundary conditions are 
?i+(0,w) = /^(WjTi) and n_{L,ui) = fB{uj,Tfi). Distri- 
bution functions relax towards the average distribution 
no = (n-|_ -|- n-)/2 with relaxation time r. Mode disper- 
sion in a one-dimensional chain is iLi{q) = 2cjo sin(ga/2), 
where a is the lattice constant. Note that the disper- 
sion of the discrete chain is used in Eqs. (|37ap and 
(j37bp as usual (see, e.g., Ref. [1^). Mode velocity is 
v{uj) = duj / dq = awo^/l ~ (w/2a;o)^ and the density of 
states is D{oj) = dq/du! = v{oj)^^. 
The solution of BTE is 



n+ix,io) = f{io,TL)-Ciij) 



2A(w) 



.ix,u)^f{u,TL)~C{oo) 



(38a) 
(38b) 



where A{uj) = tv{oj) is the scattering length and C{Ld) = 
[/b(c^,Tl) - /B(w,Tfl)]/[l + L/2A(a;)]. The solution re- 
sults in the heat current (again for h = ks = 1) 



duj 



2n 



ujv{uj)D{uj)[n+{x,uj) — n-{x,uj)] (39) 



doj 



27r l + L/2A(a;) 



[./B(w,rL) - fB{uJ,TR)], 



(40) 



which can be interpreted as Landauer-Biittiker current 
with the effective transmission function 



T 



//(^) = - 



1 



L/2K{uj)' 



(41) 



For L ^ A (a;) and classical statistics (T ^ wq), the 
thermal conductivity n = limAT-s-o derived from Eq. 

coincides with the expression obtained for the clas- 
sical self-consistent heat bath model [l^, H^l when the 
BTE relaxation time t is identified with the inverse of 
the bath coupling constant 7. This shows the similarity 
of SCTB and BTE models in infinite classical systems. 
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III. SOLVING THE SELF-CONSISTENT 
EQUATIONS 

The bath temperatures in the center region are deter- 
mined by demanding that the average heat current {Qi) 
to the local heat baths i S {1,2,..., Nc} given by Eq. 
(PU)) vanishes. Since the bath temperatures appear in 
the Bose-Einstein functions, the equations are non-linear 
and the temperatures must be solved by using iterative 
methods or by resorting to linearizing approximations. 
We use both approaches and compare solutions obtained 
from the full non-linear equations with linear and classi- 
cal approximations. Solutions of the non-linear equations 
are calculated using the Newton-Raphson method and in 
some cases an integration method based on the existence 
of a steady state towards which the system evolves. 

The Newton-Raphson method has previously been 
used to solve the SCTB equations for ID chains [23 and 
is quite an efhcient and reliable method for solving more 
general problems as well, especially when the linearized 
solution is used as the initial guess. However, each itera- 
tion of the Newton-Raphson method requires evaluating 
N'^ frequency integrals, which makes the method heavy 
for large systems. 

A slightly different and potentially better-scaling 
method for solving the equations can be found by writing 
a set of equations for the time evolution of the temper- 
atures of the local baths and letting the system evolve 
towards the steady state. If the reservoir is imagined to 
have heat capacity Ci , the temperature Ti of the reservoir 
changes due to the in-flow or out-flow of thermal current 
and obeys the differential equation 

^ = ^(Q.m,...,r^J), (42) 

where the time t is now macroscopic time such that any 
fluctuations in Qi vanish in the timescale of interest. This 
is an ordinary differential equation (ODE) of first or- 
der that evolves toward a steady-state where the bath 
temperatures satisfy the self-consistent temperature con- 
dition (Q,(Ti,...,rjVc)) = Vi e {l,2,...,iVc}. We 
call the method of integrating Eq. (|^^ the ODE method. 
In addition to having an intuitive physical interpretation 
as transient time evolution of the heat bath tempera- 
tures, ODE method has the advantage that at each time 
step, one only needs to calculate Nc frequency integrals 
to calculate the time-derivative {dTi/dt, . . . ,dT]\r^/dt). 
For large systems, the method could therefore provide a 
good alternative to the Newton-Raphson method. The 
heat capacity Ci simply affects the time-scaling in Eq. 
(|42| and can be included in the time variable t. 

The full solution of the nonlinear equations can be 
avoided by two common approximations that provide a 
linear set of equations for the bath temperatures (26j . 
Linear response approximation is based on the assump- 
tion that temperature differences are small, allowing one 




semi-infinite ]\r semi-infinite 

chain chain 

(b) 

semi-infinite lead, temperatures determined semi-infinite lead, 
temperature Tl self-consistently temperature Tr 




FIG. 2. (Color onhne) Illustration of the systems studied in 
Sees. IIV Al and IIV Bl (a) a chain of length connected to 
two semi-infinite chains, (b) a constriction of width W'^ and 
length L*^ between two leads of width W'". layers of atoms 
in the leads are included in the self-consistent calculation to 
account for the gradual temperature drop near the constric- 
tion. 

to make in Eq. ([^5]) the substitution 

fBiuJ,Tj) - /b(c^,T,) ^ ^(^,T„)(rj - T,), (43) 

where the mean temperature is Tm = {Tl + Tfl)/2. The 
linearization typically produces too low bath tempera- 
tures compared to the exact results [2^. By considering 
a single-site model, we have traced this feature back to 
the fact that the second derivative of the Bose-Einstein 
function with respect to temperature is strictly positive. 

Unlike the linear-response approximation, the classical 
approximation 

fs [io, Tj) - fs (^, T,) ^ - {Tj - Ti). (44) 

makes the self-consistent equations linear also in lead 
temperatures, so the scaling of lead temperatures by a 
constant simply scales the self-consistent bath tempera- 
tures by the same factor. 

Both linearizations exclude any non-linear effects such 
as thermal rectification j26| and produce more symmetric 
temperature profiles than the non-linear equations due to 
the equivalence of mapping Ti ^ Tl + Tr — Ti and the 
spatial reflection of the structure. 

IV. NUMERICAL RESULTS 

To highlight the pertinent physics and the properties of 
the exact and approximate solutions of the self-consistent 
equations, we study in more detail the thermal conduc- 
tion and temperature profiles in two structures shown in 
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Fig. [21 In the one-dimensional setup of Fig. [DJa), the 
temperatures are determined self-consistently in the cen- 
ter region consisting of a chain of N atoms. The chain 
is connected to two semi-infinite chains interacting with 
heat baths at constant temperatures Tl and Tr so that 
there is no geometric scattering and phonon flow is re- 
duced only by interactions with the local heat baths. The 
setup reduces to the Rubin- Greer geometry (2^ if the 
heat baths are removed. 

In the two-dimensional constriction geometry of Fig. 
[5l^b), two wide leads are connected by a narrow con- 
striction. The center region includes not only the con- 
striction but also layers of lead atoms to account 
for the effects of temperature drop near the constric- 
tion. In both geometries, nearest neighbors are assumed 
to be connected by harmonic springs with spring constant 
k ~ itiu!q, where m is the mass of the atoms. Each atom 
has only a single degree of freedom corresponding to, e.g., 
the atomic displacement in the out-of-plane direction. 

Unless otherwise stated, we set = 1 in the follow- 
ing so that dimensionless temperatures are in units of 
huJo/kB and thermal currents in units of HuJq. The di- 
mensionless friction parameter is then in units of wo- The 
friction parameter in the leads is set equal to the fric- 
tion in the central region, 7c = 7l = 7_R = 7- In Sees. 
IIV AlUVBl and UVCi all exact self-consistent temperature 
configurations are calculated using the Newton-Raphson 
iteration with the linear response temperatures used as 
the initial guess. Newton-Raphson and ODE methods 
are compared in Sec. IIVDI 



A. Rubin-Greer chain 

Due to the lack of geometric scattering in the Rubin- 
Greer setup of Fig. Elja), the setup serves as an ideal sim- 
plified model to compare the basic differences and simi- 
larities between the exact and approximate solutions of 
the self consistent problem. Figure [3] compares the self- 
consistent quantum exact, quantum linear and classical 
bath temperature profiles in a chain of length iV = 5 
with friction parameters (a) 7 = 10~^ and (b) 7 — O.I. 
The lead temperatures are set to Ti = 0.2, Tj^ = 0.1. In 
the nearly ballistic system of Fig. [3lja), all temperature 
profiles are nearly constants as a function of position, 
because coupling to baths is too weak for efficient ther- 
malization. For increased damping in Fig. [3fb), there is 
a clear temperature gradient due to interaction with the 
heat baths. The temperature gradients of the quantum 
exact and quantum linear response models are approx- 
imately the same, but the classical gradient is clearly 
larger. The most prominent feature in Figs. [3lja) and 
is, however, that the quantum exact temperature is 
higher than the temperatures obtained in linear approx- 
imations. as noted also earlier for Ohmic leads |28| . 

Figure m shows the thermal current Q = {Qr) through 
the chain as a function of left lead temperature for 
fixed right lead temperature Tr = 0.2. Friction parame- 
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FIG. 3. (Color online) Bath temperature profiles in a chain 
of length N = 5 sandwiched between two semi-infinite leads 
at temperatures Tl = 0.2 and Tr = 0.1. Friction parameters 
are (a) 7 = 10~^, i.e. the system is nearly ballistic and (b) 
7 — 0.1. Symmetry requires that T3 = 0.15 for the linearized 
models, but the quantum exact temperature is higher due to 
the non-linearity of the Bose-Einstein function. 



ters are set to (a) weak friction 7 = lO^'' and (b) strong 
friction 7 = 0.1. The length of the self-consistently mod- 
eled chain is iV = 10. In the ballistic limit of Fig. St^a), 
the current flowing through the chain at low bias ~ Tfl 
is equal to Q = Gq{Tl — Tr)^ where the quantum of ther- 
mal conductance jjl is Gq = irk'^T/Gh, which reduces 
to Gq — ttT/Q in present units. When friction is in- 
creased in Fig. m^b) , currents decrease due to the phonon 
damping caused by the heat baths. With both weak 
and strong friction, the classical approximation strongly 
overestimates thermal current, but the linear response 
approximation is valid up to ^ 0.6. The classical ap- 
proximation also makes the current response fully linear. 

Figure [5] compares the thermal currents given by the 
self-consistent thermal bath (SCTB) solution and Boltz- 
mann transport equation (BTE) solution (|40p as a func- 
tion of chain length N. In the SCTB model, the friction 
parameter in the center region is 7 = 0.1 and in the leads 
iL — iR = 0.1 or 7i = 7/f = 0.001. The string length 
L in BTE is set to L — {N -f l)a to correspond to a 
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FIG. 4. (Color online) Thermal current Q as a function of 
lead temperature Tl for fixed Tr, = 0.2. Chain length is 
N = IQ and the friction parameters are (a) 7 = 10"^ and 
(b) 7 = 0.1. With both weak and strong friction, the linear 
response approximation reproduces the exact current up to 
very high values of bias. In (a), current Q = Gq(Tl — Tr) 
corresponding to the quantum of thermal conductance Gq — 
ttT/G with T = 0.2 is also shown (black dashed). 



chain of atoms and the relaxation time in the chain is 
T = = 10. The temperatures of the baths are set to 
Tl = 0.2 and Tr = 0.1 so that the system is in the non- 
linear low-temperature regime. As expected, the current 
decreases in both models as a function of chain length due 
to phonon decay in the chain. The BTE result matches 
the exact SCTB result perfectly, when the lead friction 
parameters jl and 7^ are small, i.e., the leads are as- 
sumed to be nearly ballistic and the center region friction 
parameter is tied by the relation 7 = r^^. The require- 
ment of ballistic leads is natural, since we assumed that 
the phonon occupation in the leads is given by the Bose- 
Einstein distribution, which in SCTB model is exactly 
valid only in the limit of zero broadening, 7^ = 7l — s- 0"''. 
We have verified that the BTE and SCTB heat currents 
for small 7^ = 7^ agree also at other temperatures. Fig- 
ure [5] also shows that increasing and 77? in the leads. 
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FIG. 5. (Color online) Exact thermal current as a function of 
chain length A''. The friction parameter in the center region 
is 7 = 0.1. The BTE current (|40ll matches the self-consistent 
current if the leads are nearly ballistic, 7^=7^ = 0.001. The 
bath temperatures in the left and right semi-infinite chains are 
Tl = 0.2 and Tr — 0.1, but the currents also agree at other 
temperatures for ballistic leads. 



which increases scattering, slightly reduces the thermal 
current flowing through the center region. 

Despite the similarities between the predictions of 
BTE and SCTB for the simple ID geometry, the mod- 
els are not equivalent. For more complex geometries, 
the Green's function method, which contains full atom- 
istic dynamics, wave effects and geometric scattering, is 
a drastic improvement over solving BTE under gray ap- 
proximation. 

The agreement of thermal currents between the SCTB 
model and BTE would have been very cumbersome to 
highlight, if the leads had been described by Ohmic 
reservoirs as in earlier works instead of Rubin-Greer 
chains. Because Ohmic baths at the ends of the chain 
would reflect some of the phonons back to the chain, the 
thermal boundary conditions for the distribution func- 
tions of right and left-moving phonons in the BTE for- 
mulation would have been different from simple Bose- 
Einstein functions. 



B. Constriction in two-dimensional lattice 

In real systems, phonon transport is more complicated 
than in a one-dimensional chain due to, e.g., phonon re- 
flections from boundaries. The new features arising from 
mode mismatch at contacts and other geometric factors 
will be studied in the constriction geometry of Fig. Eljb), 
where the atoms are set in a square lattice such that a 
constriction of width and length connects two 
leads of width W'^. To account for the effects of tem- 
perature drop near the junction, atom layers in the 
leads closest to the junction are also included in the self- 
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FIG. 6. (Color online) Bath temperature profiles in a W — 
5, = 9 constriction coupled to leads of width = 71 and 
= 35 [see Fig. [IJb)]. Lead temperatures are Tl ~ 0.2 and 
Tr = 0.1. Figures show (a) quantum exact and (b) classical 
self-consistent bath temperature profiles. Friction parameter 
is 7 = 0.01. The separation of isolines is 0.05 and four contour 
lines are labeled for convenience. 



consistent calculation. The junction geometry has been 
studied earlier using molecular dynamics [i^ll^l, but in 
contrast to molecular dynamics, the present methodology 
allows to include full quantum statistics in the phonon 
populations. 

Figure [6] shows the (a) quantum exact and (b) classical 
temperature profiles in a W'~^ = 5, L'~^ = 9 constriction 
coupled to leads of size ~ 71 and = 35 in the 
low-temperature (Tl = 0.2, = 0.1) and nearly bal- 
listic regime (7 = 0.01). The asymmetry arising from 
the non-linearity of the self-consistency equations is very 
prominent in the quantum exact profile of Fig. Eta), 
as the temperature profile patterns in the left and right 
sides are visibly different. The junction temperature is 
approximately 0.17, which is notably higher than the av- 
erage temperature 0.15. This is a similar effect as noticed 
in the previous section for the ID chain: the mixing of 
the statistics of phonons at hot and cold temperatures re- 
sults in a thermal population whose temperature is higher 
than the average temperature. 



In the classical case of Fig. [HJb), on the other hand, 
symmetry of the self-consistent model requires that the 
temperature profile is symmetric with respect to spatial 
reflection and mapping Ti — >■ T^ + Tn — Ti. Therefore, the 
central part of the junction is at temperature 0.15. In ad- 
dition, the temperature profile in the bulk parts exhibits 
directional features at 45 degree angles with respect to 
the junction. These features have been observed also ear- 
lier for simil ar g eometry in classical molecular dynamics 
simulations [43|. In the quantum profile, these diagonal 
directional features are absent and the temperature pro- 
files are more directed straight towards the leads. This 
feature is even more prominent for narrower constrictions 
(not shown). The difference between the quantum and 
classical profiles is most likely related to the transmis- 
sion properties of high-energy phonons (w > T), whose 
populations are overestimated by classical statistics. An- 
other major difference between quantum and classical 
statistics is that the currents fiowing through the struc- 
ture are Q = 84.1 x 10^"* for quantum statistics and 
Q ~ 138 X 10"'^ for classical statistics, i.e., current is 
very strongly overestimated by the classical statistics in 
the low-temperature regime, as noted also for ID chain 

[Fig. m. 

Our results indicate that the diagonal temperature pat- 
terns observed in Fig. [SJb) and in the classical molecular 
dynamics simulations of Ref. (47| may be washed out 
by the quantum effects at low temperature. At higher 
temperature, quantum effects are reduced and the di- 
agonal features reappear, but only if phonon transport 
remains close to ballistic. Increasing the temperature 
also increases phonon-phonon scattering j3l|, so finding a 
temperature regime where classical statistics prevail but 
phonon transport is suflaciently ballistic can be problem- 
atic. 

The self-consistently modeled center region of Fig. [6] 
contains 4873 atoms. For this size of system and tem- 
perature range, the determination of the quantum linear 
response temperature profile, which was used as the ini- 
tial guess for Newton-Raphson iteration, took approxi- 
mately five hours wall-time with 12 CPU cores. Newton- 
Raphson iteration converged after three iterations and 
took approximately 14 hours wall-time. The calculation 
of the classical temperature profile took approximately 
14 hours wall-time as well. The solution of the classical 
temperature profile is computationally more demanding 
than calculating the linear response profile, since the pop- 
ulation functions appearing in the equations decay more 
slowly and need more integration time. 

Note that even though the system is smaller than 
the mean free path of long-wavelength phonons, the use 
of non-reflecting boundary conditions (i.e., semi- infinite 
leads) ensures that phonons are not reflected from the 
ends of back to the junction. If the ends had been 
thermalized with Ohmic heat baths, reflections from the 
baths could skew the temperature profiles. 
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FIG. 7. (Color online) (a) Graphene nanoconstriction. The 
leads extend infinitely to the left and right, but the tempera- 
tures are determined self-consistently only for the gray atoms 
in the shown center region, (b) and (c) Self-consistent bath 
temperature profiles (K), in (b) quantum exact case and (c) 
classical approximation. The semi-infinite leads are at tem- 
peratures Tl ~ 300 K and Tr = 280 K. The relaxation time 
r = 1/7 is set to 1 ps. 



C. Thermal transport in a graphene constriction 

The two-dimensional square lattice model is easily ex- 
tended to real materials such as graphene. It is inter- 
esting to see, for example, if the directional features ob- 
served in the square lattice remain for more complex lat- 
tice geometries. The example geometry is shown in Fig. 
[Tlja). Each atom now has three degrees of freedom, which 
are all coupled to a single local Langevin heat bath. We 
set the temperature range close to the room tempera- 
ture, Tl = 300 K and — 280 K, because the acoustic 
phonon lifetime r at room temperature is known to be 
of the order of t = 1 ps suggesting that the bath 
coupling constant is 7 = = 10^^ s~^. Carbon-carbon 
interactions are modeled by the fourth-nearest-neighbor 
force constant model [i^ with the parameters of Ref. 



[Hoj . which reproduce the bulk ab-initio phonon spec- 
trum of graphene very accurately, at least for the acoustic 
modes. The optical modes are not active at room tem- 
perature, since they are populated only at temperatures 
close to T « 1000 K, but are fully included in the model 
in any case. We have verified the correct implementation 
of the force constant model by comparing ballistic ther- 
mal conductances of pure nanoribbons to the results of 
Ref. dH. 

Figures [TJ^b) and [TJ^c) show the quantum and clas- 
sical temperature profiles close to the room tempera- 
ture. The temperature profiles agree quite closely, which 
is unexpected since the phonon populations originat- 
ing from the classical and Bose-Einstein distributions at 
room temperature are quite different: the highest-lying 
vibrational energies of graphene correspond to temper- 
atures of Td « 2300 K. The agreement of temperature 
profiles therefore suggests that only the low-frequency 
modes close to the F point, for which quantum and clas- 
sical statistics agree, contribute to the transport and de- 
tailed temperature profile. On the other hand, the heat 
flow through the structure is still quite strongly overesti- 
mated by the classical approximation: quantum current 
is Q « 2.1 X 10-8 W and the classical Q « 5.1 x 10"^ W. 
No directional features appear in the studied geometry at 
room temperature, but lowering the temperature and in- 
creasing the phonon relaxation time could produce more 
complex temperature profiles. These studies are left for 
future work. Note also that approximately only half of 
the total temperature drop takes place in the constric- 
tion. 



D. Comparison of the solution methods 

As a final example of the numerics of the solutions, 
we compare the Newton-Raphson iteration and the or- 
dinary differential equation (ODE) method. The differ- 
ential equation ([42|) is integrated using the MATLAB® 
(5^ implementation of an explicit Runge-Kutta formula 
with the Dormand-Prince pair (ssj and adaptive step- 
size. Using an adaptive step-size integrator is neces- 
sary to avoid slowing down as integration approaches the 
self-consistent temperature configuration. Integration is 
stopped when the maximum heat current flowing to the 
bath is less than 10-^7. 

Figure [5] shows the comparison of the two methods in 
the search for self-consistent solution. The setup is the 
Rubin-Greer setup of Fig. [2ta) with N = 2 and the val- 
ues of the friction parameter are (a) 7 = 0.01, and (b) 
7=1. Although it is best to use, e.g., the linear response 
approximation temperatures as initial guess for iteration, 
we use now Ti = Tl = 0.2 and T2 = Tr = 0.1 for illus- 
trative purposes. The contour lines of the target function 
f{Ti,T2) = Ql + Ql are also shown. The function / is 
defined such that the self-consistent temperature config- 
uration is the global minimum and zero of /. The self- 
consistent temperatures are {Ti,T2) = (0.1590,0.1585) 
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FIG. 8. (Color online) The search for the self-consistent 
bath temperatures Ti and T2 that satisfy Qi(ri,r2) — 
Q2{Ti,T2) = for an TV = 2 chain. The lead tempera- 
tures are Tl ~ 0.2 and Tr =0.1 and the friction parame- 
ters are (a) 7 — 10~^ and (b) j — 1. The two methods used 
for the search are Runge-Kutta-Fehlberg integration of Eq. 
(|42|) (blue dots connected with dashes) and Newton-Raphson 
iteration (red crosses connected with dash-dots). The con- 
tours belong to the "target function" /(Ti, T2) that is defined 
to be the squared sum of the currents flowing to the self- 
consistent reservoirs, /(Ti.Ta) = Qi{Ti,T2f +Q2{Ti,T2f. 
The self-consistent temperature configuration {Ti,T2) satis- 
fies f(Ti,T2) — and is also the global minimum of /. The 
initial guess is Ti = 0.2, T2 = 0.1. 



and {fi,f2) = (0.1616,0.1574) for the cases of Figs. E^a.) 
and[8jb), respectively. 

For both -weak and strong friction, the Newton- 
Raphson iteration proceeds similarly: The first iteration 
step of the Newton-Raphson method slightly misses the 
solution, but the second iteration already takes temper- 
atures very close to the self-consistent temperature con- 
figuration. ODE method, on the other hand, proceeds 
approximately along the direction of steepest descent in 
the target function. Since the gradient V/ = 2J-^Q, J 
being the Jacobian matrix of Q, is not necessarily paral- 



lel to Q, the path taken by the ODE method is generally 
not strictly along the steepest descent. 

For the case of larger 7 in Fig. ^h), the contour lines 
are elongated forming a canyon-like shape and the heat 
exchange between the local baths starts to dominate over 
the heat exchange with the leads. In this case, ODE 
method does not proceed directly towards the solution. 
We have noted that such cases can be very difficult to 
handle for the ODE method, since the residual time inte- 
gration along the canyon requires a very small step size. 
Newton-Raphson iteration, on the other hand, always 
seems to find the solution with only a few iterations. 

In future, it would be interesting to study how well 
phonon damping and dephasing induced by the self- 
consistent heat baths mimics true anharmonic effects. 
Comparisons could be carried out, for example, by com- 
paring the classical approximation of self-consistent equa- 
tions with classical molecular dynamics (MD) simula- 
tions. Knowing now that the bath temperatures cor- 
respond to the local kinetic temperature [Eq. 
the local bath temperatures could be meaningfully com- 
pared to the local kinetic energy densities obtained from 
MD. It would also be worth investigating whether non- 
diffusive transport effects such as anomalous heat con- 
duction in one dimension (s^ could be reproduced by 
using a frequency-dependent bath coupling constant. If 
that is the case, one could study quantum effects in 
anomalous transport using the Green's function method. 



V. CONCLUSIONS 

We studied quantum heat transport in nanostructures 
using the Green's function method combined with self- 
consistent heat baths. Semi-infinite leads acting as ther- 
mal reservoirs were reduced to sources of noise and 
dissipation in the boundaries of the scattering region. 
In the scattering region, the temperatures of the heat 
baths mimicking anharmonic effects were determined 
self-consistently from the requirement of heat current 
conservation. The self-consistent bath temperatures were 
shown to measure local energy density, thereby giving 
them a meaningful physical interpretation. In the classi- 
cal limit, local kinetic temperature is equal to the bath 
temperature. 

By coupling one-dimensional chain to semi-infinite 
chains, thereby eliminating contact resistance, we demon- 
strated the equivalence of thermal currents obtained by 
the self-consistent thermal bath model and the Boltz- 
mann transport equation under gray approximation with 
full phonon dispersion. Self-consistent thermal bath 
model is, therefore, a physically meaningful method to 
introduce phonon relaxation to ballistic quantum trans- 
port models. 

As an application of the formalism, we presented tem- 
perature profiles in two-dimensional constrictions and 
showed that quantum statistics plays a vital role in how 
directional patterns of temperature emanate from the 
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junction. In a graphene constriction at room temper- 
ature, the bath temperature profile obtained by classi- 
cal approximation agreed very closely with the quantum 
temperature profile, suggesting that quantum effects are 
not strong and molecular dynamics simulations could be 
justified under those assumptions. In more general cases, 
we expect, however, that the Green's function method 
combined with self-consistent thermal baths is a very use- 
ful tool in studying quantum heat transfer in the ballistic, 
diffusive and crossover regimes of phonon transport due 
to the good balance of complexity, insight and predictiv- 
ity it offers. 
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Appendix: Derivation of Eq. I|29p 

In this appendix, we derive Eq. ([29| for thermal cur- 
rent flowing to a local heat bath. For notational sim- 



plicity, we combine the lead noises 77^ and ijn and center 
region local bath noises ^c* to a single vector variable 
s-^, where the index J £ {L, R,l,2, . . . , Nq} is now a 
general index for either a lead bath (J G {L, R}) or a 
self-consistent local bath (J G {1,2, Nc})- Nq is 
the number of atoms in the center region. Explicitly, 
= ?/L , = rjn and s* is a vector whose only non-zero 
component is s!- — ^cij the ith component of ^c- The 
noise covariances are then 

{s'' {uj)s''' {uj')) = 2TrS{oj + oj')r' (uj) [fBiL^,Tj) + 1]S'^'^' . 

(A.l) 

If index J G {L,R}, the coupling function F"' is defined 
by the self-energy of the lead, Eq. ([^5]) . For J = i, the 
only non-zero matrix element in the coupling function is 
F*j(aj) = 277170; (we write 7(7 = 7 in this section). 

Equation (fT3|l can be written (dropping the index C 
for center region) 

u{uj) = -G{uj) s^ioj). (A.2) 
J 

The heat flowing to an Ohmic bath is obtained by cal- 
culating the statistical average of the symmetrized heat 
current 

= 777iu^ - ^[ui£,i + Ciiii]. (A. 3) 

We proceed term by term. The statistical average of the 
first term is 



i^mu',) = 7m Q ^i-iLu)ud^)e-''^' J ^(-jw')u.(c^')e"*"'*) (A.4) 

= / ^^7m(-c.a;')e-''"+"''*EEG'.^-(^)G^'»(^')(^i(^)^^''(^')> 
J ^ jj, jk 

= I ^7"^u;2E [fB{i^,Tj) + 1] . (A.6) 
The average of the second term is 

-(^.6) = - (/ '^{-^u:)u.{u:)er^-' j ^^.l^Oe"^"'*) (A.7) 

= / S^^-*")^'^^"''"'^*EE^-(") (^/(-)e.(-')) . (A.8) 

J j 

The only term surviving the sum over baths J is the one corresponding to the local heat bath at site i, so 

/I 

^iGuiiohmuj^ [fB{oo,T,) + 1] . (A.9) 

Combined with the symmetrizing term, one gets 

-7^{i^,i,+^^u^) = - / -^AGuiuj) ~ G,, (-w)]77nc^2 [/^(^^T,) + 1] (A.IO) 

= -J ^E[<^Hr''HG(-^)]n7rnu;2[/s(^,T,) + l], (A.ll) 



where we used Eq. with the replacements g/ ^ G and f ^ X)/ ■ Combming Eqs. (|A.6|) and (jA.ll[) . we get 
(Q.) = I I^T^^'E [C^(^)r'(^)G(-c^)1,, [Ib{u^.Tj) fB{i^,m . (A.12) 
Noting that the integrand is an even function finally gives Eq. ([29]) . 



[1] S. Berber, Y.-K. Kwon, and D. Tomanek, Phys. Rev. [29 

Lett. 84, 4613 (2000). 
[2] P. Kim, L. Shi, A. Majumdar, and P. L. McEuen, Phys. [30 

Rev. Lett. 87, 215502 (2001). 
[3[ C. W. Chang, D. Okawa, H. Garcia, A. Majumdar, and [31 

A. Zettl, Phys. Rev. Lett. 101, 075903 (2008). 
[4[ L. G. C. Rego and G. Kirczenow, Phys. Rev. Lett. 81, 

232 (1998). [32 
[5[ K. Schwab, E. A. Henriksen, J. M. Worlock, and M. L. 

Roukes, Nature 404, 974 (2000). [33 
[6[ M. Prunnila and J. Mehaus, Phys. Rev. Lett. 105, 

125501 (2010). [34 
[7[ I. Altfeder, A. A. Voevodin, and A. K. Roy, Phys. Rev. [35 

Lett. 105, 166101 (2010). [36 
[8[ A. Majumdar, Science 303, 777 (2004). 
[9[ Y. Dubi and M. Di Ventra, Rev. Mod. Phys. 83, 131 [37 

(2011). 

[10] E. Pop, Nano Res. 3, 147 (2010). [38 
[11[ M. Terraneo, M. Peyrard, and G. Casati, Phys. Rev. 

Lett. 88, 094302 (2002). [39 
[12[ B. Li, L. Wang, and G. Casati, Appl. Phys. Lett. 88, 

143501 (2006). [40 
[13[ C. W. Chang, D. Okawa, A. Majumdar, and A. Zettl, 

Science 314, 1121 (2006). [41 
[14[ J. Y. Murthy, S. V. J. Narumanchi, J. A. Pascual- 

Gutierrez, T. Wang, C. Ni, and S. R. Mathur, Int. J. [42 

Muhiscale Comp. Eng. 3, 5 (2005). 
[15| R. Landauer, Philos. Mag. 21, 863 (1970). [43 
[16| M. Biittiker, Phys. Rev. B 46, 12485 (1992). [44 
[17| N. Mingo and L. Yang, Phys. Rev. B 68, 245406 (2003). [45 
[18| N. Mingo, Phys. Rev. B 74, 125402 (2006). 
[19| J.-S. Wang, J. Wang, and N. Zeng, Phys. Rev. B 74, [46 

033408 (2006). [47 
[20[ J.-S. Wang, J.Wang, and J. Lii, Eur. Phys. J. B 62, 381 

(2008). [48 
[21[ M. Bolsterh, M. Rich, and W. M. Visscher, Phys. Rev. 

A 1, 1086 (1970). [49 
[22| F. Bonetto, J. L. Lebowitz, and J. Lukkarinen, J. Stat. 

Phys. 116, 783 (2004). 
[23[ W. M. Visscher and M. Rich, Phys. Rev. A 12, 675 [50 

(1975). 

[24[ A. Dhar and D. Roy, J. Stat. Phys. 125, 801 (2006). [51 
[25| D. Roy, Phys. Rev. E 77, 062102 (2008). 

[26| D. Segal, Phys. Rev. E 79, 012103 (2009). [52 
[27[ E. Pereira, H. C. F. Lemos, and R. R. Avila, Phys. Rev. 

E 84, 061135 (2011). [53 
[28] M. Bandyopadhyay and D. Segal, Phys. Rev. E 84, 

011151 (2011). [54 



R. J. Rubin and W. L. Greer, J. Math. Phys. 12, 1686 
(1971). 

U. Weiss, Quantum Dissipative Systems, 3rd ed. (World 
Scientific, Singapore, 2008). 

J. Ziman, Electrons and Phonons: The Theory of Trans- 
port Phenomena in Solids (Oxford University Press, 
USA, 2001). 

G. W. Ford, J. T. Lewis, and R. F. O'ConneU, Phys. 
Rev. A 37, 4419 (1988). 

M. P. Lopez Sancho, J. M. Lopez Sancho, J. M. L. San- 
cho, and J. Rubio, J. Phys. F: Met. Phys. 15, 851 (1985). 
N. Li and B. Li, J. Phys. Soc. Jpn. 78, 044001 (2009). 
J.-S. Wang, Phys. Rev. Lett. 99, 160601 (2007). 

G. W. Ford, M. Kac, and P. Mazur, Journal of Mathe- 
matical Physics 6, 504 (1965). 

A. Dhar and B. S. Shastry, Phys. Rev. B 67, 195405 

(2003) . 

C. Caroli, R. Combescot, P. Nozieres, and D. Saint- 
James, J. Phys. C: Solid State Phys. 4, 916 (1971). 
T. Yamamoto and K. Watanabe, Phys. Rev. Lett. 96, 
255503 (2006). 

W. Zhang, T. S. Fisher, and N. Mingo, Numer. Heat 
Transfer, Part B 51, 333 (2007). 

H. Goldstein, Classical Mechanics (Addison- Wesley Pub- 
lishing Company, Inc., 1971). 

G. D. Mahan, Many Particle Physics, 3rd ed. (Kluwer 

Academic/Plenum Publishers, New York, 2010). 

S. Lepri, R. Livi, and A. Politi, Phys. Rep. 377, 1 (2003). 

A. Majumdar, J. Heat Transfer 115, 7 (1993). 

A. J. Minnich, G. Chen, S. Mansoor, and B. S. Yilbas, 

Phys. Rev. B 84, 235207 (2011). 

S. K. Saha and L. Shi, J. Appl. Phys. 101, 074304 (2007). 
K. Saaskilahti, J. Oksanen, R. P. Linna, and J. Tulkki, 
Phys. Rev. E 86, 031107 (2012). 

N. Bonini, J. Garg, and N. Marzari, Nano Lett. 12, 2673 
(2012). 

R. Saito, G. Dresselhaus, and M. S. Dresselhaus, Physical 
Properties Of Carbon Nanotubes (Imperial College Press, 
1998). 

L. Wirtz and A. Rubio, Solid State Comm. 131, 141 

(2004) . 

Z. Huang, T. S. Fisher, and J. Y. Murthy, J. Appl. Phys. 
108, 094319 (2010). 

MATLAB, 8.0.0.783 (R2012b) (The MathWorks Inc., 
Natick, Massachusetts, 2012). 

J. Dormand and P. Prince, J. Comp. Appl. Math. 6, 19 
(1980). 

A. Dhar, Adv. Phys. 57, 457 (2008). 



